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in different one-dimensional geometries, whereby tunneling, resonance and anharmonic- 
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Chebyshev expansion technique of the time evolution operator. 
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Introduction 

Quantum statistical physics, such as condensed matter or plasma physics, but also 
quantum chemistry, heavily depends on effective numerical methods for solving 
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complex few- and many particle problems. Implementing suitable theoretical con- 
cepts for their description on modern (super-) computer architectures, nowadays 
computational physics constitutes, besides experiment and theory, the third col- 
umn of contemporary physics ^ . Numerical techniques become especially important 
for strongly correlated systems where analytical approaches largely fail. This may 
be due to the absence of small (coupling) parameters or, more general, because the 
relevant energy scales are not well separated, both preventing the application of 
standard perturbative schemes. Another challenging problem, that calls for numer- 
ical approaches, concerns the description of microscopic and nanoscopic systems, 
and of particles in finite quantum structures (restricted geometries), where the level 
quantization might become as important as particle correlations. 

A first idea for the numerical study of this kind of quantum systems might be a 
brute force exact diagonalization of the underlying (model) Hamiltonian. Consider- 
ing the exponential growth of the Hilbert space with the number of particles, such 
a description of quantum many-particle systems is (and will be in the future) out 
of reach. The situation becomes even more difficult if bosonic degrees of freedom 
(e.g. phonons in a solid) come into play, resulting in an infinite dimensional Hilbert 
space even for finite systems. One way to circumvent this problem is to restrict the 
many-particle Hilbert space to the physically most important subset. Along this 
line, e.g., density matrix renormalization group schemes have been developed EIH 
Semiclassical descriptions offer another possibility to overcome these limitations. A 
variety of semiclassical methods has been proposed during the last decades, espe- 
cially to describe the dynamics of quantum systems 4 5 6 7 8 9 xhese methods are 
appealing, in some sense, as their close relation to concepts familiar from classical 
physics facilitates an intuitive interpretation of the results. Also in view of bridging 
the gap between quantum many-particle and classical continuum theories they seem 
promising for describing systems in the thermodynamic limit. Therefore we will fo- 
cus on semiclassical approaches to quantum dynamics in the following, of which we 
will discuss three carefully selected ones in more detail. 

The majority of the traditional semiclassical methods is based on a real time 
path integral formulation of quantum mechanics [MIH Expressing the time evolu- 
tion of the complex wave function in terms of action integrals, the Feynman path 
integral is generally evaluated within the stationary-phase approximation. Then the 
occurring integrals can sometimes be performed analytically, otherwise numerically 
using direct integration or Monte Carlo (MC) techniques Integrating an oscil- 
latory, complex valued integrand, the dynamical sign problem, however, spoils to 
some extend the efficiency of the MC integration. Despite the exponential decay 
of the integrand outside a vicinity of the classical trajectories the numerics is still 
challenging. In the sum of many contributions with different complex phases most 
of them may cancel out and the result may become exponentially small. 

Equivalently, a quantum system can be described in terms of real valued quan- 
tum phase space distribution functions'^, e.g., the Wigner function This over- 
comes the problem of handling a complex valued wave function but, because of 
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possible negative values of the Wigner function and the Heisenberg uncertainty re- 
lation, the Wigner function cannot be interpreted as a joint probability. Instead, 
it should be considered as a convenient mathematical tool for the description of 
quantum systems. In the numerical work, the dynamical sign problem is alleviated 
for the Wigner function but still present l^^ l ^^ l i'^i 

In order to overcome the heavily debated dynamical sign problem, some years 
ago the description of quantum states in terms of a positive function, the so called 
quantum tomogram, has been proposed The strict positivity of the tomogram 
seems promising in view of an effective MC sampling of the trajectories. In the 
framework of the tomographic representation a description of quantum dynamics 
by diffusive Markov processes was suggested in Recently, the applicability of the 
tomographic approach to one- and two-particle systems has been demonstrated 
However, to the best of our knowledge, most of those studies are based on harmonic 
potentials but validations for arbitrarily shaped potentials are still missing. 

Motivated by this situation, it is the aim of the present work, to outline and com- 
pare the various methods quoted above with respect to their accuracy and compu- 
tational performance. The selected semiclassical methods were chosen to represent 
algorithms based on different descriptions of quantum states: the wave function, 
the Wigner function and the quantum tomogram. Instead of investigating complex 
physical systems in which a multitude of effects compete, we benchmark the meth- 
ods on the basis of relatively simple systems, concentrating onto a single aspect in 
each case. We will focus on the dynamics of a wave packet in distinct model geome- 
tries which account for basic aspects of quantum mechanics: tunneling, confinement 
and nonlinearity (anharmonicity). Calibrating the different approaches by studying 
simple toy models is necessary in order to detect their limitations and prospects, 
before applying them to the more complicated problems of current interest, e.g., 
to classical chaotic systems, quantum localisation, dynamic tunneling or entangle- 
ment. The obtained semiclassical results will be compared with an exact solution 
for the time evolution of the quantum system. Thereby, the exact solution is calcu- 
lated using a highly efficient expansion of the time evolution operation in a series 
of Chebyshev polynomials. 



1. Computational schemes 

1.1. Chebyshev expansion of the time evolution operator 

As a reference for the approximate results we will present in the next sections, 
the Chebyshev expansion of the time evolution operator constitutes a very efficient 
technique which fully includes all quantum effects. Governed by the time dependent 
Schrodinger equation, the dynamics of a quantum state |'0(to)) in a time indepen- 
dent external potential may be expressed in terms of the time evolution operator 
U{t,to) = U{At) = e-''H{t-to)/h ^-^^Yi At = t - to. Expanding U{At) int o a serie s of 
first kind Chebyshev polynomials of order /c, Tk{x) = cos(A: arccos(x)), ESEMI, we 
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obtain 

M 

U{At) = e-^^^^/^ " 



co{aAt/h) + 2 V Ck{aAt/h)Tk{H) 



k=l 



(1) 



Prior to the expansion, the Hamiltonian has to be shifted and rescaled such that 
the spectrum of H = {H — b)/a is within the definition interval of the Chebyshev 
polynomials, [—1, 1]. The parameters a and b are calculated from the extreme eigen- 
values of H ash = ^(£^max + ^min) and a = ^(^max — ^min + e). Here, we introduced 
e = a(£^max — Emin) to cusurc the rescaled spectrum |^| < 1/(1 + a) to be well 
inside [—1,1]. In practice, we use a = 0.01. Note, that the Chebyshev expansion 
also applies to systems with unbounded spectra. In those cases we truncate the 
infinite Hilbert space to a finite dimension by restricting the model on a discrete 
space grid or using an energy cutoff. By this, we ensure possibly large but finite 
extreme eigenvalues. 

In ([T]), the expansion coefficients Ck{aAt/h) are given by 

Ck{aAt/h) = / \ ^ dx = {-ifJk{aAt/h) , (2) 

J TTV 1 — 



where Jk denotes the k-th order Bessel function of the first kind. 

To calculate the evolution of a state |'^(to)) from one time grid point to the 
next one, = /7(At) |?/;(to)), we have to accumulate the c^-weighted vectors 

\vk) = Tf^{H)\tlj{to)). Since the coefficients Ck{aAt/h) depend on the time step but 
not on time explicitly, we need to calculate them only once. The vectors \vk) can then 
be calculated iteratively using the recurrence relation of the Chebyshev polynomials 

l^fc+i) = 2H\vk) - \vk-i) , (3) 

where \vi) = H\vo) and \vo) = |V^(to)). 

In the numerics, we use a discrete coordinate space basis \qi), i = 1, . . . , TV with 
{Q\Qi) = ^{Q — Qi)^ representing an equally spaced grid of points. The wave function 
at position qi is given by the i-th entry of the corresponding (complex) vector, 
i^iQi^to) = {Qi\i^{to))' Aiming at the description of a spatially unbounded system, we 
choose the extension of our space grid such that ?/;(gi, t) = ilj{qN^t) = throughout 
the whole simulation. In this way, no artificial reflections at the boundaries of the 
simulation volume arise and the results are independent of the actual size of the 
simulation volume. Overall, evolving the wave function from one time step to the 
next requires M matrix vector multiplications (MVM) of a given complex vector 
with the (sparse) Hamiltonian of dimension N as well as the summation of the 
resulting vectors after appropriate scaling: 

M 

m)) = e-''""'/^ Jo{aAt/h)^2^{-i)'MaAt/h)n{H) |V^(to)) . (4) 

k=l 

Note that the Chebyshev expansion may also be applied to systems with time 
dependent Hamiltonians. However, there the time variation H(t) determines the 
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maximum At by which the system may be propagated in one time step. For time 
independent in principle, arbitrary large time steps are possible at the expense 
of increasing M. 

1.2. Linearized semiclassical propagator method 

Aiming at the description of a many-particle system, the numerical effort for a full 
quantum calculation increases drastically. If the system is described in terms of 
(correctly symmetrized) product states, the main numerical problem is the expo- 
nential growth of the Hilbert space dimension with the number of particles. For 
some systems, the full correlations encoded in those states are of minor importance 
as particular aspects of interest may be described on a lower level of complexity. 
Therefore, much interest has been devoted to finding a suitable semiclassical ap- 
proximation for the time evolution operator. Based on a path integral description of 
quantum mechanics 113111, the Suzuki- Trotter decomposition HHHI of U{t — to) opens 
the road toward a class of semiclassical approximations. Instead of directly propa- 
gating the quantum state as a whole using U{t — to)^ we consider the propagation 
of several individual paths from qo{to) to q{t) (virtual particle trajectories). The 
corresponding propagator is given by 

U{q,t]qo,to) Ceyip{iS{q,t]qo,to)/h) , (5) 

V. paths 
qorvq 

with some normalization factor C and the action 

t 

S{q,t;qo,to) = J dt' p{t')0) - H W),q(t')] , (6) 

^0 

evaluated along each virtual path. H{p^ q) is the classical Hamilton function of the 
system. The only restriction for the virtual paths are fixed starting and end points, 
go(^o) and q{t)^ apart from which they are completely arbitrary. Specifically, they 
do not follow the Hamilton equations of motion, which only hold for a subset of 
them, namely the classical paths (left panel of Fig. [T]). 

In a numerical implementation (middle panel of Fig. [T]), the quantum state is 
represented on a discrete coordinate grid qi for any time grid point, while the coordi- 
nates of the virtual trajectories are treated as continuous variables. Reconstructing 
the quantum state at time t involves the integration over the propagators between 
all possible combinations of q^ and q. While we choose the go from the set of dis- 
crete grid points, the end points q may also be off the grid. For the deposition!^ 
of a trajectory to the qi grid, a suitable shape function A{qi^q), e.g., a Gaussian or 
simply a delta peak at the nearest grid point, is necessary, 

i^{Qi^t) = J dqo J dq A{qi,q)U{q,t;qo,to)il^{qo,to) . (7) 

As the contributions to ip^qi^t) from the different paths are complex valued, their 
superposition includes interference effects in the reconstructed quantum state. The 
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Fig. 1. Left panel: Virtual (blue, dashed-dotted lines) and classical (red, solid line) paths con- 
necting qo(to) and q(t). For the classical path the action along the trajectory is extremal. Middle 
panel: Visualization of the deposition scheme ([7|. At each coordinate grid point qi the initial wave 
function ijj^qi^to) determines the sampling weight for all trajectories starting from qi. The path 
propagation follows the classical equations of motion and the virtual (non-extremal) paths are 
taken into account perturbatively. As the end points of the paths q{t) need not match the grid, 
they are deposited there using the shape function A(q,qi). Right panel: Wave function propaga- 
tion as implemented in this work. Around each grid point the potential is approximated to first 
order. The propagator for a constant force field is known analytically. Starting from all possible 
pairs of initial conditions (qi,Pj), i, j = 1, . . . N, the trajectories are deposited in momentum space, 
matched to the momentum grid by A(p,pj). The reconstruction of the wave function in coordinate 
space is done by inverse discrete Fourier transform. 



major difference between existing semiclassical propagator methods concerns the 
numerical implementation of The standard line of argumentation in the lit- 
erature is that the main contributions to (|5| are those for which S is extremal - 
these are just the classically realized trajectories - and the paths in their vicinity. 
This requirement is deduced from stationary phase integration. Then, the action is 
expanded to second order around the extremal trajectories giving the well known 
WKBFI^ result M 



cpsiths ' 
qorxq 

where Sdq^t] qo^to) is the action along the extremal trajectories from qo{to) to q{t) 
and the phase of is fixed to 7r/4. The sum accounts for the fact that specifying 
go and q does not uniquely determine a classical trajectory. In addition to the 
phase factor from the classical action, the phase of the propagator is determined 
by the sign of d'^ Sd {dq^dq). Using po = —dSc/dqo, the Morse theory IMIMI allows 
for separating the determinant of the monodromy matrix dq/dpo from its sign. 
Relating the number of focal points (at which dq/dpo vanishes) along a trajectory 
to the number of negative eigenvalues of d'^Sc/ (dqdqo), the whole phase information 
can be encoded in the Morse index u. For each focal point, u is increased by one, 
and the square root is taken from the absolute value of dq/dpo. Then we obtain the 
Van Vleck Gutzwiller propagator EHHl 



c. paths 
q^r^q 



dq 



1/2 

^iSc{q,t;qo,to)/h^-mu/2 ^ j^gj 



dpQ 

In view of a numerical evaluation, this formulation has several shortcomings. First, 
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the so called root search problem consists in finding all initial momenta po for which 
trajectories starting in go end up in q. For complex systems, this boundary value 
problem is much more demanding than the solution of an initial value problem. 
Second, at the focal points the expression (|9| diverges and the expansion of S 
has to be taken to higher order. At those points it is particularly difficult to keep 
track of the correct branch of the square root of dq/dpo and thus determining 
the Morse index. Both problems are circumvented in contemporary applications 
of the semiclassical propagator by expressing ([9| in an initial value representation 
(SCIVR) i^i^ i ^^ i Inserting ^ into ([7|, this is achieved by a change of variables, 

j dqo J dq j dgo j dpo 



c. paths 



dpo 



(10) 



circumventing the root search, and the monodromy matrix appears in the numera- 
tor, 

1/2 



= -L= I dqo I dpo A{quq{qo,po)) 



dq{qo^Po) 



dpo 



(11) 



Here, the final coordinate q{qo^Po) is a deterministic function of the initial values qo^ 
Pq. Instead of using a coordinate space basis, similar methods have been formulated 
in momentum or coherent state representation IMIllES, leading to the Herman-Kluk 
propagator EIEH For all these methods the problem of the oscillatory integrand 
still persists. To improve the convergence properties of the MC integrals, various 
integral filtering techniques have been proposed in the literature EE^EIIIH Thereby, 
the basic idea is to filter out the high frequency oscillations of the integrand which 
contribute only little to the integral but are the main obstacle for an efficient MC 
evaluation. 

In this work, we follow a recently proposed slightly different route (cf. 

right panel of Fig.jl]). Discretizing the potential on the coordinate grid g^, we locally 
approximate the potential to first order and use this approximation for qi — Aq/2 < 
q < qi -\- Aq/2. A justified question then is, to which extend the concatenated 
potential will be able to fully describe quantum effects and reproduce the results 
for the original continuous potential. The propagator for a trajectory in a linear 
potential is known analytically 



m 



27rih{t - to) 
X exp 



h \ 2{t-to) 2' ' 24m 



(12) 



where s is the slope of the potential. In contrast to the WKB formulation, here 
tracking the sign of the monodromy matrix is not necessary, but the maximum 



8 



possible time step is severely reduced. This limit is set by the validity of the local 
potential approximation, i.e., by the potential variation and grid spacing. The time 
step has to be chosen such that even the fastest particles stay within their initial 
grid cells intermediately. An advantage of the constant force field approximation is 
the bijection between final positions and initial momenta which allows for a straight 
forward formulation of the algorithm in terms of an initial value representation. The 
range of initial momenta, which are required by Fourier completeness is given by 
the coordinate grid spacing as \po\ < irh/Aq. A final trick consists in depositing 
the contributions from all virtual particles not directly to the coordinate grid. In- 
stead, they are gathered in momentum space and the reconstructed wave function 
in coordinate space is obtained using inverse discrete Fourier transform, 

v-i*'*) = E -p (^) ^fe'*) ■ (13) 

j = l V TT \ / 

This mixed representation has proven less noisy than a direct deposition in coordi- 
nate spacel^. In analogy to ([7|), we use the shape function A{pj^p) for the deposition 
on the momentum grid 

ij{pj,t) = J dqo J dp A{pj,p)U^'''{p,t;qo,to)ij{qo,to) , (14) 



where go, ^o) is the Fourier transform of (12), 



n''^{p,t;qo,to) = -i= / d9^""(<^,^;<^o,^o)e-*^''/^ (is) 



Instead of an explicit evaluation of the g-integral in ([15|, we profit from the knowl- 
edge about the trajectories, 

q{t) = -^{t-tof + ^{t-to) + qo, (16) 
Zm m 

to change the integration variable from q to po. Finally, this allows us to propagate 
the wave function by one time step. At = t — to, using an initial value representation. 



X exp 



i f m[q{qo,po) - qo]'^ sAt a , i s^{At)^ 



h\ 2At 2 L^^^"'^"^ 24m 



b(<7o,Po) + ^o] — pq{qo^Po) 



(17) 



where q{qo^Po) is given in (16). 



1.3. Wigner-Moyal approach 

Dating back to Wigner the idea of describing quantum mechanics without a 
complicated operator algebra, but by equations for commuting variables, has been 
very appealing. Then, introducing a phase space distribution function and replacing 
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operator expressions according to the corresponding rules of association, quantum 
expectation values may be calculated by simple integration over commuting vari- 
ables. One of the most common distribution functions is the Wigner function, 

W{q,p, t) = ^J e-'^''r{q " Vn/2, t)^P{q + i^h/2, t) dr? , (18) 

together with the corresponding Weyl rule of association Within this rule of 
association, operator expressions ordered as 6^^*^+^^^ are replaced by their scalar 
variables, e.g., e^^^+^^^ ^ ^i^g+ir^p^ with ^, 77 G C. Starting from the von Neumann 
equation for the time evolution of the density matrix, we determine the evolution 
equation for the Wigner function l^^l^'^H^I as 



dt m dq 



dW f 

i^(^)^= J dsW{q,p-s,t)co{s,q). 



(19) 



where F{q) = — is the classical force, and 

„,,,) = _||d<,'n,-,>in(?f:)+i^«,)!f^. (20, 

In the classical limit, the right hand side of (19) vanishes, leaving us with the 
Liouville equation for the phase space density. Then the dynamics can be expressed 
in terms of the classical propagator 

n'^((7,p,t;go,Po,^o) = 8[q - q{t]pQ,qQM)\^[p - Pit'.Po.qoM)] • (21) 

Here p and q are the momentum and coordinate of a trajectory which evolves 
according to the Hamilton equations of motion subject to the initial conditions 
p(^o) = Po and q{to) = qo. Using H^, we may rewrite (19) in form of an integral 
equation mmEU 

W{q,p,t) = J dpodqo {q,p,t; qo,Po,to)Wo{qo,po,to) 

t CO 

J dr J dprdqr {q,p,t; qr,Pr,r) J ds W{qr,Pr - s,t)u;{s, qr) , 

to —CO 

(22) 

and solve it by iteration. To lowest order, we only keep the first line in ( [22| ) and 
neglect the second integral completely. This means, we propagate classical trajecto- 
ries {q^p) in time, after sampling their initial conditions po and qo from the initial 
Wigner function Wo(^o,Po,^o) at time to using a MC procedure. Assembling those 
contributions at the next time grid point t = to + At, we obtain the lowest order 
approximation W^^\q^p^t). 
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Fig. 2. Cartoon of propagating one trajectory in presence of momentum jumps. At time to, the 
initial conditions for the trajectory (go,^>o) are sampled from VKo(go, Po, ^o) by a MC procedure. 
Then up to time r G [to,^] this trajectory is propagated to {qr^Pr — s) according to the classi- 
cal equations of motion. There, an instantaneous change in momentum by s occurs, leaving the 
intermediate position qr fixed. From the new phase space point {qr^Pr) the trajectory evolves 
again classically up to following the propagator {q^p^t^qr^PT^T). When summing up the 
contributions of all trajectories, each one has to be weighted by the function C(;(s,gr), accounting 
for the momentum jump s and the potential at the intermediate position qr. 



For the second order approximation, t), we expand the Wigner func- 

tion in the last integral in ([22|) consistently to lowest order, 



f w (2^) 

= / dpodqoU {Qr.Pr - s,r;qo,Po,to)Wo{qo,Po,to) . 

This allows for evaluating ly^^^ also by means of trajectory methods. A sketch 
of the basic idea behind the second order approximation is given in Fig. [2j together 
with a detailed description in the caption. Including momentum jumps in the evo- 
lution of W^^) we take into account that continuous phase space trajectories are 
guaranteed only for classical or almost classical systems Depending on the mo- 
mentum jump s and the potential V{qr — q')^ the weighting factor uj{s^qr) may 
become negative. Therefore, sign changes of the Wigner function at some phase 
space points are included in the second order approximation, whereas they are ab- 
sent in the first order approximation due to strictly positive trajectory weights. 



Using the improved estimate W^'^^ in (23), the construction of higher order terms 



is straight forward. Conceptually, higher order terms correspond to the inclusion 
of several momentum jumps within one time step. As the fundamental aspects 
(discontinuous trajectories, possibility of negative trajectory weights) are already 
included in W^'^\ we restrict ourselves to the two lowest orders of the approximation. 
A detailed investigation of the inclusion of higher order terms is beyond the scope 
of this work (for details see J^H^. 
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Having W{q^p^ t) at hand, it is not difficult to calculate expectation values of all 
kinds of operators from it. For any operator expression which is ordered as e^^^+^^^, 
we follow the Weyl rule of association and replace each operator by its corresponding 
scalar function. After multiplication with the Wigner function we integrate over the 
whole phase space and obtain the desired expectation value. 

1.4. Tomographic approach 

In view of a probabilistic interpretation of quantum mechanics, the Wigner function 
has the shortcoming of possible negative values, which prevents its interpretation as 
a probability density. Furthermore, any interpretation of a joint probability, that si- 
multaneously determines coordinate and momentum for a quantum system, violates 
the Heisenberg uncertainty relation and is therefore misleading. Meaningful results 
for probabilities and expectation values, are integrals over the Wigner function in 
extended phase space regions, e.g., minimum uncertainty Gaussians, or along phase 
space contours. Such a contour integration is used in the tomographic representa- 
tion of quantum mechanics proposed some years ago The so called quantum 
tomogram Ell^ 

w{X,fi,u,t)= I ^^^^^W{q,p,t)e-'^^^-^^-^P\ (24) 
J 27r 

relates to the Wigner function by a class of Radon transformations'^ which are char- 
acterized by II and v. A simple, intuitive interpretation of the quantum tomogram 
is shown in Fig. [sj^a). The vector (/i, v) fixes a direction, with X the distance from 
the origin. Then the tomogram w{X^ ja^ v) is just the integral over the Wigner func- 
tion along a straight line perpendicular to (/i, v) which passes through X. Choosing 
(/i, v) appropriately, we may continuously change between coordinate and momen- 
tum representation. Describing a quantum system in terms of the tomogram is 
completely equivalent to the wave function or Wigner function representation. Its 
relation to the Wigner function is obvious due to its definition, and we obtain the 
Wigner function from the tomogram by the inverse map 

W{q,v,t) = I ^^^^(X,^,^,t)e^(^--''--) . (25) 

The relation to the wave function formalism is a little more involved. We may extract 
the probability densities in any rotated reference frame (in particular coordinate 
and momentum space) from the tomogram, but not the wave function itself. The 
phase information inherent to the wave function is distributed over all reference 
frames (/i, u). Each fixed choice of (/i, v) gives only a density information, e.g., the 
coordinate representation \ip{q)\'^ is equal to w{X, ^ = 1, = 0). Keeping this aspect 
in mind as well as the larger required set of variables, one might ask at this point 
what the profit of this method shall be. Considering the dynamics of a quantum 
system, the strict posit ivity of the tomogram is its main numerical advantage. This 



12 





P 




(n(to),v(to)) 














W(q,p,to) 







P 

(H(t),v(t)) 










W(q,p,to) 


"<W(q,P,t) 



(^(t'),v(t')) 




W(q,p,t 



Fig. 3. Cartoon of the quantum tomogram and its time evolution. Panel (a): The tomogram 
it;(X(fo), m(^o), ^(^o),^o) is the integral over VK(g,p,to) along the thick red line. Panel (b): Extrac- 
tion of the tomogram in the same reference frame at time t by evolving the Wigner function and 
keeping (/ir, ^r) = [l^it)^ ^(^)] = [m(^o), ^(^o)] fixed. Panel (c): Time evolution of the tomogram by 
the method of characteristics. Keeping the Wigner function W{q,p^ to) fixed, (/ir, i^r) = [l^(t), i^{t)] 
is propagated backward in time to {/~i{t'), ^(^0) using \29\ for different local approximations. Sum- 
ming all those contributions gives the tomogram at t in the desired reference frame (/ir, i^r)- 



circumvents the dynamical sign problem encountered in the Wigner function and 
semiclassical propagator approaches. 

The time evolution of the tomogram can be thought of in two equivalent ways. 
For any reference frame (/i(to), ^^(to))^ e.g., the coordinate representation (1,0), the 
tomogram can be calculated for all times by keeping (/i, u) fixed and evolving the 
Wigner function in time [Fig. [3]^b)]. This is exactly the interpretation we used in 
Sect. 1.3 to extract t)p from the Wigner function. Alternatively, we can exploit 



the evolution equation of the tomogram .22i^ 



dw 
~dt' 



ji dw 
m dv 



V 



d 



dii d/dX 



\hv d 



\hv d 



(26) 

which can be derived from the von Neumann equation. In view of the definition 
of w in (24), the 'anti-derivative' d^^w just multiplies W{q^p) by i/k. Therefore, 
the term d^d^w is well defined. A solution of (26) for harmonic potentials V{q) = 
^mujQ{q — Qc)'^ can be given explicitly. In analogy to the continuity equation for the 
tomogram. 



dw dw dw 



dt 



dt dX 



dw . dw . 



(27) 



we collect the terms in (26) and find the correspondences 



X = mu?{ 



oQc 



v = —ii/m. 



(28) 



This set of linear differential equations defines trajectories in (X, /i, v) space charac- 
terizing the reference frames in Fig.js] Solving the system (28 ), we get the propagator 
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from (Xo,/io,^^o,^o) to as 



1^0 



sin[cJo(^-^o)] - i^ocos[cjo(^-^o)] 



x8 



X - Xq- /io<?c(cos[cc;o(t-to)] - l) - i^o^^o<7c sin[cc;o(t-to)] 
ji — jiQ cos[a;o(^ — ^o)] + vomcoo sin[a;o(^ — ^o)] 



(29) 

The uniqueness of the z^(t)) trajectories is due to the g'-independency of 



the coefficients in the set of equations (28). 



Evaluating (26) for arbitrary potentials (beyond the free particle or harmonic 
case), higher order derivatives of w appear. This spoils the identification with the 



continuity equation (27) and the application of the method of characteristics (28). 



A possible way out is the local expansion of the potential to second order. Then 



(26)-(29) are valid locally for each q with appropriate parameters uJo{q) and Qciq) 



which are determined from the harmonic expansion around q. Since now several 



propagators (29) exist, the [X^ji^v) trajectories are not unique anymore, and the 
tomogram at the new time grid point is given by 



w{X,ii,v,t) = J dXo J d/io J d^o J dg 



(30) 



^loiq),qciq) ^' ^' -^0, MO, i^O, ^o)^(^0, MO, ^0, ^o) • 



Usually, we are not interested in the full w{X^ (i^v^t) but only in a particular 
reference frame {/j^r^i^r)^ e.g., the coordinate representation. Having the method of 
characteristics in mind, we search for those trajectories (X(to), m(^o), ^(^o)) that 
give contributions to {X(t)^ ii(t)^v{t)) = {Xr^ /J^r^^r)- To this end, we start from 
the known iD(X^, /i^, i^^, to) and evolve the trajectories backward in time to f = 
to — {t — to) according to (29) using different coordinates q [cf. Fig.jsj^c)]. The desired 
w{Xr^ /J^r^^r^ t) is then the sum over all such tomograms w{Xq(t')^ l~^q{t')i ^^(^0? ^o), 
where the index q reflects the dependency on the parameters uoQ^q) and qc{q) in the 
propagator 

An alternative approach which relates the time evolution of the quantum to- 
mogram to the theory of Markov processes has been considered recently 
Interpreting the propagator 11^ as a transition probability for a Markov random 
process, its time evolution is governed by the Chapman-Kolmogorov equation*^ 



n^(z,t;zo,to) = j dzr U^{z,t;Zrjr)U^{zr,r;zo,to) , 



(31) 



where to < r < t and we introduced the abbreviation z = {z^^z'^^z^) := (X, /i, z^). 
Due to the positivity of the tomogram and its normalization, the requirements 



U^{z,t;zo,to) >0 , J dzoU^{z,t;zo,to) 



1 



(32) 
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for a Markov process are fulfilled. Furthermore, the locality in time of the Hamil- 
tonian guarantees that no memory effects are present. The dynamics of diffusive 
Markov processes may be equivalently described by two partial differential equa- 
tions, the first and the second Kolmogorov equation 

3 ™t 3 3 



dU7 

dto ^ 

dt 



0, 



i=i j^i 

3 3 



dz'dz^ 

i=l j = l 



in which the drift and diffusion coefficients are defined as 

a'{zr,r)= lim — ^ — / dzr' {zr' Zr,r){zl, - z^) . 
t'^t r - t J 



b^^{zr,r) = lim 



dzr' U^{zr',r]Zr,r){zl, - zl){z^^, -4)- 



(33) 
(34) 

(35) 
(36) 



In the limit r' ^ r we can evaluate ( 35 ) and ( 36 ) using the harmonic approximation 
n^J^ from (29). As the harmonic propagator depends on uJo{q) and qdq) of the 
local potential expansion, also the drift and diffusion terms are q dependent. For 
clarity of the notation, we will suppress the additional index for the moment. 

Instead of solving (34) directly, we consider the equivalent system of stochastic 
integral equations [MIl3for the underlying random variables z. 



z\t) = z 



3 

to to 



(37) 



Here ^-^(r) are white noise random processes with (^-^(r)) = and (^-^ (T)^^(r')) = 
S{r — t')5^^ . For evaluating the second integral we will refer to the Stratonovich 
definition of a stochastic integral The functions ilj{zr^r) and g'^^{zr^r) relate to 
the drift and diffusion coefficients in (35) and (36) as 



a\zr,T) = tp'izr.r) + 2 

j=l k^l 

3 

b^^{Zr,T) = J29'''{Zr,T)g^''{Zr,T). 
k=l 



dZr 



(38) 



(39) 



For an implementation of the derivative term in (38) we profit from the equiva- 
lence!^ 

j = l k^l / 

where (...) means the stochastic expectation value over the normalized Wiener 
process ^/e(r) and A^/e(r) is the corresponding increment of this random process 
in Ar. As a consequence of the Stratonovich integration the matrix elements g'^^ 
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have to be evaluated at 



Zr + ^Azr, i.e., shifted by half the increment of Zr 
during Ar. According to (39) these matrix elements can be calculated via Cholesky 
decomposition 1^ from b{zr,r). Note, that the matrix b is only positive semidefinite, 
requiring some care in the numerical determination of ^, e.g, adding e tob and taking 
the limit e ^ 0. 

In summary, the evolution of a (X, /i, u) trajectory is calculated iteratively from 
(37)- (40). To lowest order, in ( [37| ) only the deterministic drift term (35) is taken 
into account, which is evaluated using the local harmonic propagator (29). From 
the resulting increments z{t) — z{to) we get a first estimate for b (and thus g) which 
we use in the next iteration. Having access to the deterministic part as well as the 
diffusion term in (37) we calculate several trajectories, needed for the expectation 
value in (40). Now an approximation for all terms in (37) is available and we can 
finally propagate our (X, /i, z/) trajectory. Alternatively, we can repeat the last step 
to ensure that the increments entering in (40) have been calculated using the full 
stochastic integral equation. 

In addition to the stochastic character of the trajectory propagation, we have to 
keep in mind that the drift and diffusion terms depend on the coordinate at which 
(29) is evaluated. A suitable importance sampling of q is of crucial importance for 
an efficient implementation. 

Although we discussed for simplicity only a three-dimensional vector z = 
(X, /i, u)^ which corresponds to one set of initial conditions, the generalization to k 
initial conditions directly carries over. In any case, including the initial condition 
for the coordinate representation is obligatory as a reconstruction of q out of z is 
necessary for intermediate evaluations of (29). 

From the tomogram, general expectation values can be calculated in analogy to 
the Wigner function case If the desired expectation value involves only position 
or momentum operators but not both, this task simplifies to an integration over the 
corresponding density. Let us consider for instance the kinetic energy ^{p^)- With 
jii = and f = 1^ the integral representation of the ^-distribution in (24) shows that 
X = p and we get 



1 

2m 



^/dXX^^(X,, = 0,. = l) = ^/ 



(41) 



2. Numerical Evaluation 

As test cases for the above methods, we consider one-dimensional systems described 
by the Hamilton operators Hi = ^ ^Vi (q) . The four potentials 

ViiQ) = ^'^^W + Vbexp(-g^) , V2{q) = ]:mujlq^ - Vbexp(-g^) 

^1 ^1 (42) 

= 2^^o(<7^ + (^sq^) , V^iq) = Vq ^ -mujl{-q^ + a^q"^) , 

shown in Fig. [4j each pose another numerical difficulty to semiclassical approaches. 
It is well known, that semiclassical methods perform well for harmonic or weakly 
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anharmonic cases, but that the effects of strong anharmonicities are hard to cap- 
ture Specificahy, we chose Vi{q) to study tunnehng effects, V2{q) in view of 

resonances and Vs{q) to investigate the influence of a nonhnear force term. Finahy 
V4(g) combines the chahenge of tunnehng effects and anharmonicities in the poten- 
tial. The inclusion of the shallow harmonic trap in Vi{q) and V2{q) prevents the 
particle from escaping to infinity after the scattering event and restricts the simu- 
lation volume to a reasonable size. For all cases we use the same initial conditions, 
a Gaussian wave packet of width a, centered at go, with center momentum po, 

^{q.t = 0) = ^^J^y/A {"^(^ " ^o)' + ^^^4 * ^^^^ 

Taking the modulus squared \^l^{q^t = 0)p of ( [431) , correct normalization to 
unity is obvious. Using the initial wave function ( |43[ ) together with (18) we get the 
corresponding initial Wigner function 

M^(g,p,t = 0) = ^exp|-^(g-go)'-^(p-Po)'| ■ (44) 

Comparing the coefficient of {p — Po)'^ in the exponent to the standard form of a 
Gaussian, the standard deviation in momentum space reads ap = h/ {2a). Therefore, 
aap = /i/2, which makes our initial state a minimum uncertainty Gaussian wave 
packet. Integrating the Wigner function over the momentum (coordinate) variables, 
we get the probability density in coordinate (momentum) space, 

mQ,t = 0)\^ = I dp W{q, p,t = 0) = -jL= exp | - | . (45) 

It is clear, that we could have obtained this expression also directly by taking the 
modulus square of (43). To obtain the probability in momentum space, 

|^(p,t = 0)|2 = |dgH^(q,p,t = 0) = y^exp|-?^^fa^|, (46) 

from the wave function, however, first a Fourier transform to momentum represen- 
tation is necessary. 



i,{p, t = Q) = -i= / Aq e-'f V(9, t = 0) 
V27ra J 



2a^V'^ \ a\ ,2 i i 



(47) 



(48) 



Finally, we obtain the tomograms of the initial state using (24) as 

v/27rcr2(/i,i/) I 2cr^(/i,i/) J 

where the width of the tomogram gt depends on the particular reference frame 
(/i, v\ and is given by 




ar{ii,v) = \\v'^\ — \ +m2<t2. (49) 
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Fig. 4. From left to right: Benchmark potentials Viicj) as given by ( |42| . The initial state in all 
cases is the same. For the first case we show the real and imaginary part of the wave function by 
dashed and dashed-dotted lines. In the other panels, the modulus square of the wave function is 
given. For comparison of the relevant energies, the baseline of the wave function is drawn at the 
energy of the initial state. 



Throughout this work we express ah quantities in terms of the fixed reference 
units for length (ui)^ mass {um) and time (ut). From those we may also construct 
units for energy {ue = Umuf/u^) and momentum {up = UmUi/ut). In these units 

Um, u;oUt = 0.1, asuj 0.01, 
As initial conditions, we choose qo = —5u£, 



h = Umuj/ut. To be specific, in (42) we use m 
uo/^Ut = 0.4, a/^u\ = 0.02 and Vq 
Po — Up and a = ugl\f2. 



2.1. Discussion of the time evolution 

Probability densities Figure |5] shows the time evolution of the probability density 
in coordinate space, for times up to t = 56ut. Each column corresponds 



to one of the above potentials (42) while the rows refer to the methods used. The 



topmost row, calculated by Chebyshev expansion (C), gives the exact solution. For 
the barrier case [Vi{q)^ first column] the particle hits the barrier (center position 
marked by a dashed line), at about t ^ 5ut. Here, the main part of the wave packet 
is refiected, but a sizeable fraction also penetrates through the barrier. The quan- 
tum nature of the particle gives rise to the interference pattern on the left hand 
side of the barrier, where high and low probability densities alternate. This effect 
is still more pronounced around t ^ 35ut when the transmitted and refiected parts 
interfere after their reunion. For the case of a quantum well [^2(^)7 second column], 
we observe the overall picture expected for the dynamics in a simple harmonic trap. 
Since the width of the initial Gaussian, however, does not match the resonance fre- 
quency of the harmonic trap, its width changes during the time evolution, refocusing 
once each tt/ujo. The rather small effect of the additional dip at g = consists in 
a reduced density in this region for all times and a retardation of the transmission. 
For the case of an anharmonic potential [Vs{q)^ third column], besides the above 
mentioned broadening of the wave packet, the nonlinearity of the force causes ad- 
ditional interference effects. Elongations up to q ^ 15ii£, exceeding the classical 
turning points qc ^ 8.2u£, are possible due to the quantum nature of the initial 
state. Despite its fixed total energy, the initial wave packet contains contributions 
with higher and lower energies. The forth column shows the rich structure of the 
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Fig. 5. From left to right: Time evolution of the probability density in coordinate space, \ip{q, 
for a particle in the benchmark potentials Vi{q) to V4(q) as given by ( |42| . For each case, we show 
(from top to bottom) the exact solution, calculated by Chebyshev expansion (C), the results from 
the linearized semiclassical propagator method (P), the results from the first order Wigner-Moyal 
approach (W) and the tomographic results (T). For implementation details see Sect. |2.2] 



dynamics in the double well potential. While the principle part of the wave function 
remains in the left well, a considerable amount penetrates the barrier where strong 
interference effects arise. The back scattering of these contributions causes strong 
interference patterns also in the left well for t > lOut. 

In the second row, the results from the linearized semiclassical propagator 
method (P) fully coincide with the exact results. The excellent agreement con- 
firms, that within a semiclassical framework it is in principle possible to capture 
quantum effects, although by an tremendous increase of computational resources 



(cf. Sect.jM). 

Restricting the iteration series for the Wigner function to first order, we obtain 
the results shown in the third row (W). The initial splitting of the wave packet for 
the barrier case Vi{q) at about t ^ 5ut is reproduced while the second one around 
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t ^ 35ut is missed. The explanation of this failure is simple. Trajectories with high 
enough energy to cross the barrier once will be able to cross it any time, oscillating 
forth and back in the trap potential. Others which already failed the first passage due 
to an insufficient energy are reflected every time they hit the barrier. Hence they stay 
on the left side of the barrier for all times and in the second half period, t > 35iit, 
there is almost no spectral weight to the right of the barrier. Nevertheless some, 
exponentially rare, trajectories can be found there. Large negative initial momenta 
make some trajectories energetic enough to overcome the barrier, but lead to a phase 
shift as compared to the trajectories with positive initial momentum. Furthermore, 
the first order Wigner results fail to resolve the fine structure of the interference 
effects for Vs{q) and V4(g). In this sense, the quantum information contained in this 
approximation is limited even though in the initial Wigner function Wo{q,p^to) all 
quantum effects (to arbitrary high orders of h) are included. As we will see later, 
taking into account the second order term of the iteration series will alleviate these 
problems. 

The bottom line shows the results from the tomographic approach (T). Here, the 



positions for the potential evaluation in (29 ) are calculated from classical trajectories 



sampled from the initial distribution. While also here, the dynamics of the system 
is described qualitatively, there is nevertheless a large discrepancy to the exact 
data. Four weak points should be stressed. First, any sharp feature in is 
washed out even more than in the Wigner approach. Second, for Vi{q) and V2{q)^ 
the signatures at t ~ 5ut erroneously extend to too large negative g-values. The 
reason for this are diverging trajectories, caused by the negative curvature of the 
barrier potential. Third, the turning points in the anharmonic potentials {Vs{q) 
and V4{q)) are highly overestimated, which is an effect of the coordinate sampling 
for the potential evaluation. Forth, the tunneling in the double well potential is not 
accounted for correctly. This is due to the large range of q for which the curvature of 
the potential is negative. Hence, the amount of rejected trajectories is overestimated. 

Expectation values The semiclassical propagation of the Wigner function and 
the tomogram is intented to give an adequate description of complex many particle 
systems. Thereby, average values take the center stage. 

For the potentials Kri,2,4}(^) the initial wave packet splits in two parts. To 
account for this in terms of expectation values, we define reduced average quantities 
for the positive and negative half- axis, 

oo 

{q}- = ^ J q\i>{q,t)\Mq , {q)+ = ^ J q\i>{q,t)\Mq , (50) 

-OO 

where N± is the partial norm on the corresponding half-axis, 

oo 

N_ = J \i^{q,t)\Mq , N+ = Jmq,t)\^dq. (51) 
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Fig. 6. Expectation values for the benchmark potentials Vi{q) (left column) and V2(q) (right 
column). The upper panels show the part of the wave function norm on the negative half-axis N—. 
In the lower panels, the mean positions {q)± on the corresponding half-axis are given. We show 
the exact results from the Chebyshev expansion (C - solid lines) as well as the first and second 
order approximations from the Wigner-Moyal approach (VK^-*^) - dotted lines; Vl^(^) - dashed- 
dotted lines) and tomographic results (T - dashed lines). As the wave functions obtained by the 
semiclassical propagator method agree within numerical accuracy with the exact results from C, 
no separate curves are shown. 



For the first two test cases, we show in Fig. |6] the time evolution of the initial 
state in terms of {q)± and The norm on the positive half axis is not shown 
separately as A/'+ = 1 — For Vi(^), the constant value of N- after the ini- 
tial tunneling event {t ~ but) indicates an independent evolution of the two wave 
packets on their corresponding half- axis. Only at the refocusing point {t ^ 35ut) 
weight is temporarily redistributed between them due to interference effects. From 
{q)±, the oscillation between the barrier and the confining trap of each wave packet 
can be identified. The mean energy of the transmitted wave packet is larger than 
for its reflected counterpart as the mean position reaches larger values. Quantum 
dynamics, however, comprehends more than a simple energy discrimination of the 
constituents of the wave packet by the barrier. 

This becomes obvious when considering the first order Wigner result, Con- 
taining the classical energy discrimination only, the finite value of for t > 35ut 
is missed as already discussed in Fig. [5] The ratio of to N- within this ap- 
proximation is solely determined by the initial energy of each simulated classical 
trajectory. Those constituents of the initial state with E > Vq will overcome the 
barrier while others will not. For an estimate, we calculate the energy of the ini- 
tial state, neglecting for simplicity the influence of the barrier. Then N- (A^+) is 
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the integral over the initial Wigner function inside (outside) an ellipse defined by 
p^/{2m) + ^mujQq^ = Vq. This result agrees well with the numerical data. The 
strong underestimate of Nj^ for t > 3but within I^^^^ explains the deviation of (q)-^ 
in this range. 

Taking into account the second order term of the iteration series, W^'^\ slightly 
improves the accuracies of both the norm and the reduced expectation values; the 
relative weight of the split wave packet {N-) after t > 35ut still deviates about 
15% from the exact value. This deviation is related to the finite grid resolution, 
smearing out the exact positions during the recurrent coordinate sampling and grid 
deposition. Increasing the grid resolution, this effect can be reduced systematically. 

In the tomographic approach, the focus on expectation values does not resolve 
the problems encountered for the probability density. Around t ^ 5ut, when mainly 
the barrier region is sampled, we observe two consequences of the negative potential 
curvature. On the one hand, the fraction of trajectories which overcome the barrier 
is overestimated (too small value of N-). On the other hand, a considerable amount 
of the reflected trajectories diverges and thus the value of of is too negative. 
Apart from this the expectation values reproduce the results of C qualitatively but 
not quantitatively. Here the especially remarkable features are the larger elongation 
of (q)- and the absence of a pronounced modulation of {q)-^ for t > 35ut. 

For V2{q) (right column) the norm indicates the nearly perfect transmission of 
the wave packet since for each time there is considerable weight on one half-axis 
only. Therefore, the average positions show an almost perfect sinusoidal oscillation, 
if we consider (g)^ for t < 35ut and (q)- afterwards (thick parts of the solid lines). 
On the respective other half-axis the average values should be taken with care due 
to the low weight. The rather good agreement between the expectation values from 
C and VK^^^ in this case is not surprising. It is known from the literature that for 
harmonic potentials the exact Wigner function can be obtained by classical propa- 
gation of trajectories'^. As the dip around g = is only a moderate perturbation, 
W^^^ describes the dynamics still well. For this case, also the method of Wigner 
trajectories I ^^ I ^Q I ^^ I is applicable, in which the higher order terms of the iteration 
series are included perturbatively into a pseudopotential V{q). The trajectories then 
follow the classical equations of motion with respect to V{q) instead of V{q). For the 
construction of such a pseudopotential an approximate Wigner function is required. 
This limits a practical application of the method of Wigner trajectories to slightly 
perturbed harmonic potentials or nearly free motions Furthermore, the existence 
of V{q) is not guaranteed due to the perturbative character of this scheme. Since for 
this potential the W'^'^^ results for {q)± almost perfectly match C for the majority 
branch, the profit in considering W^'^^ is limited. Merely in the time evolution of 
N_ an improvement from VK^^^ to VK^^^ is noticeable. The tomographic results once 
more show a qualitative agreement with C. In addition to the deviations for the 
minority branch, the diverging trajectories, arising when the wave packet crosses 
the dip, perturb the results around t/uf ~ 5,35. 
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Fig. 7. Time evolution of the Wigner function for potential Vi(g). Upper panels: Exact result 
using Chebyshev expansion (C). Middle panels: Propagation of the initial Wigner function using 
classical trajectories only [first order approximation, VF^-*^^]. Lower panels: Inclusion of momentum 
jumps [second order approximation, VK^^^]. For each method snapshots at times tjut = 4, 8, 20, 36 
are shown. For comparison, the inset in the upper left panel shows the extension of a minimum 
uncertainty Gaussian. 



Wigner function Focusing on the barrier case, Vi{q)^ we compare in Fig. [7| the 
exact Wigner function, calculated from the (Chebyshev propagated) wave function 
together with definition (18), to the first and second order approximation. The ex- 
tension of a minimum uncertainty Gaussian (inset in upper left panel in Fig. [7| is 
the lower limit where the concept of a joint probability holds due to the Heisenberg 
uncertainty relation. Nevertheless, comparing W{q^p) for different approximations 
can be taken as a good quality check since any observable of the system is obtained 
as an integral over the Wigner function. Starting from the initial Wigner function 
(Gaussian, centered at g = — 5ii^, p = Up), the overall time evolution is dominated 
by a clockwise rotation of the phase space points. When the fastest contributions 
of the wave packet encounter the barrier {t = 4iit), for the first time significant 
negative values of W{q,p) occur. The initial Gaussian shape breaks up into a tri- 
angle, refiecting the low-energy trajectories which are held back by the barrier and 
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the fast contributions overcoming the barrier. The quantum nature is reflected in 
a weak interference pattern of smah positive and negative values around the main 
structure. Evolving further in time {t = Sut)^ the regions with negative weights get 
more pronounced. At t = 20ut there is a strong interference pattern of large positive 
and negative values of W{q^p) in between the two major positive portions left and 
right of the barrier. When the transmitted part returns to the left side of the barrier 
{t = 36ut) the structure remains divided with strong interference in between the 
two positive bulks. 

Considering the first order approximation, the most pronounced difference com- 
pared to the full quantum result is the absence of regions with negative values in 
W'''^\ This is clear as the initial state W{q^p^t = 0) is strictly positive and during 
the classical propagation of the trajectories their weight is unchanged. Therefore, at 
any time W'^'^^ is a superposition of positive contributions. Despite the simplicity of 
this approximation, all regions with large positive weights are in essence captured 
correctly. Regions of nearby positive and negative values in the exact solution are 
marked within W*^^^ by vanishing, or strongly reduced values (cf. q/ui G [0,10], 
p/up G [—1,1] at t = 20ut). The integral over p in this region vanishes for both 
VK^^^ and the exact solution, explaining why physically measurable quantities like 
t)P agree well despite the differences in W{q,p). 

Including the second order term of the iteration series, the presence of negative 
weights is restored in W^'^\ Even though, arbitrarily large momentum jumps are 
allowed for the trajectories, finite amplitudes of W^'^'^ are restricted to the center 
region. Outside this region, the fast oscillations of u;{s^qr) guarantee the complete 
cancellation present in the exact results. Taking the finite simulation grid into ac- 
count, the extent of this cancellation will delicately depend on the used resolution 
and accuracy of the (MC) integration. A signature of a non-perfect cancellation can 
be seen for t = 20ut at values of q/u£ G [5, 10] and large p as a series of positive and 
negative stripes. Accumulating such numerical fluctuations, the true region with 
finite amplitudes is overestimated for large times {t = 36ut). 

2.2. Details of implementation 

Chebyshev expansion In view of further applicability of the proposed methods, 
let us focus on the computational requirements in what follows. In this regard, the 
description of quantum dynamics by Chebyshev expansion of the time evolution 
operator is an extremely powerful method. Its extraordinary performance is not 
restricted to the one-dimensional test cases considered here, but has been demon- 
strated successfully for more complex, higher dimensional cases IHIIH Despite the 
exponential growth of the Hilbert space with the number of particles, also an ap- 
plication to many-particle systems, e.g., the polaron problem is within reach 
As compared to the standard Crank-Nicholson algorithm the Chebyshev expan- 
sion has two main advantages: speedup and larger accessible system sizes. First, for 
the Chebyshev expansion only MVM are required, while for the Crank-Nicholson 
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scheme, 



1 + ^iHAt/hj mo + At)) = (^1 - ^iHAt/hj mo)) , (52) 



a linear equation system needs to be solved in each time step. While for the one- 
dimensional case considered here, the tridiagonal structure of the coefficient matrix 
speeds up the calculation, in general the solution of this system is the most time 
consuming step. Speeding up the calculation by an initial inversion of the time 
independent coefficient matrix and performing successive MVM afterwards is only 
feasible for moderate Hilbert space dimensions. 

Second, the Crank-Nicholson algorithm is accurate to order (At)^, whereas the 
accuracy of the Chebyshev expansion is determined by the expansion order M. 
We may choose M such that for k > M the modulus of all expansion coefficients 
\ck{aAt/h)\ ~ Jk{aAt/h) is smaller than a desired accuracy cutoff. This is facilitated 
by the fast asymptotic decay of the Bessel functions, 

J^{aAt/h)^^=(^y for A:^oc. (53) 
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Then, for large M, the Chebyshev expansion can be considered as quasi-exact and 
thus permits a considerably larger time step. For example, using a simulation grid of 
N = 1024 sites with grid spacing Aq = O.OSu£ the necessary scaling parameters for 
the four test cases are a = 160, 161, 226, 2307, and for A: > M = 108, 108, 140, 1028 
all \ck\ < 10"-^^. These cutoffs ensure that the wave function is exact for the used 
time step At = OAut for all times. Here, 'exact' means that within numerical 
accuracy the wave function agrees with the time dependent wave function obtained 
by a full diagonalization of the Hamiltonian, 

N 

mt)) = E e-'^"*/'^|n)(n|V(t = 0)) , (54) 

n=l 

where \n) are the (time independent) eigenstates of the system and the corre- 
sponding eigenenergies. 

Besides the high accuracy of the method, the linear scaling of computation time 
with both time step and Hilbert space dimension are promising in view of further 
applications to more complex systems. Almost all computation time is spent in 
MVMs, which can be efficiently parallelized, allowing for a good speedup on parallel 
computers. 

Linearized semiclassical propagator method The considered implementation 
of the linearized semiclassical propagator method is not intended for a fast de- 
termination of an approximate solution. In this respect, more efficient flavors of 
semiclassical propagator methods can be found in the literature, where higher order 
potential terms are taken into account and the propagated trajectories are chosen 
by some kind of importance sampling. Instead, we focus on the question how close 



25 



a semiclassical approximation can be to the exact quantum solution if we let all 
concerns about computational requirements aside. To achieve the desired accuracy, 
two aspects are of prominent importance. First, we cannot choose the trajectories 
that have to be propagated solely according to the current wave function amplitude 
at the grid points. For the overall interference effects, also trajectories starting at 
grid points with low amplitudes are of importance. Already discarding trajectories 
with initial weight \^lj{qo^ to)\ < 10~^ influences the interference pattern and leads to 
a noticeable deviation from the exact results in phase and also magnitude. Second, 
the used time step has to fulfill the Courant criterion During a single time step 
any trajectory has to stay within its initial grid cell to ensure the validity of the 
local potential approximation. Fortunately, the maximum distance a trajectory may 
cover in one time step can be calculated exactly to optimize the time step. From ([l6| 
we read the displacement in one time step At as Aq = — s(At)^/(2m) -\-poAt/m. 
Substituting the largest value for \po\ = 7rh/{Aq) and requiring \Aq\ < Ag/2, the 
maximum time step is given by 



As for the reconstruction of the wave function contributions from the whole grid 
are necessary, the grid points with the largest slope s will be most restrictive for 
the time step. For Vi{q) to V4{q) the used time steps are At/ut = 5 x 10~^,5 x 
10-^,5 X 10-"^, 10-"^ on grids with N = 1024, 1024,512,256 sites and grid spacing 
= 0.125u£. Using those parameters, the results reproduce within numerical 
accuracy the results from the full quantum calculation, including the correct phase 
of the wave function. 

Wigner-Moyal-approach The numerical demands of directly propagating the 
Wigner function depend drastically on the order of the iteration series taken into 
account. For the classical propagation and assembly of a large number of 

paths {Np ^ 10^) is possible at very low computational costs. The continuity of 
the phase space trajectories allows for once sampling the initial conditions and then 
following those paths up to arbitrary times. In contrast, for higher order approxima- 
tions a single initial sampling is not sufficient anymore. Due to momentum jumps 
also regions far away from the classical end points acquire a finite weight. Then, a 
resampling of the initial conditions at each time grid point is necessary. Perform- 
ing such a resampling also for W^^'' slightly influences the data shown in Fig. ^ 
While for short times the agreement is almost perfect, with increasing time the 
sharp features present in the continuous results are washed out by the resampling. 
Caused by the grid discretization this effect systematically decreases with the used 
number of grid points. The shape function used in the deposition to the grid influ- 
ences the necessary number of trajectories to achieve a desired accuracy. Using the 
'cloud-in-ceir (CIC) scheme weight is attributed to the two nearest grid points 
in each direction, constituting a reasonable compromise between deposition costs 




(55) 
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and necessary broadening. For the first order results shown in Figs. [5]- [7| we used 
Np = 10^ initially sampled trajectories. Those are deposited by the CIC scheme 
onto a 400 x 400 grid with Aq = 0.22bui and Ap = OMbUp. 

For W^'^\ the computational requirements increase drastically. While the recipe 
how to implement this order approximation is straight forward, improving the ac- 
curacy as compared to the first order results is numerically challenging. Especially 
the stability of the long time evolution depends drastically on several factors. First, 
in each time step the MC sampling of the initial conditions depends on the pre- 
vious result. Therefore, numerical and statistical fluctuations may amplify in runs 
with too poor statistics. As for the considered one-dimensional systems a direct 
integration is feasible, we refrain from the MC integration to circumvent possible 
convergence problems. Adapting Fig. |2]to the full integration scheme, the consid- 
ered initial phase space points (go,Po) cover the whole grid. As long as t — to 
is not too large, the r-integral can be evaluated by the midpoint rule. In absence 
of momentum jumps the classical trajectory evolves continuously in [to^t]. As com- 
pared to the dependency on the magnitude of the momentum jump occurring at r 
the influence of the exact jumping time r in [to,t] is of minor importance. Work- 
ing on a fine s grid, at r all those momentum jumps are performed for which the 
final position {Qt^Pt) does not exceed the simulation grid. From the updated phase 
space point, the corresponding trajectories are then evolved up to time t, where 
they are deposited onto the (q^p) grid by means of the CIC scheme. Second, for 



the calculation of uo{s^qr) the derivative of the (^-distribution in (20) needs to be 
implemented numerically. Approximating the ^-distribution by a Gaussian of width 
cr^, the corresponding derivative reads 



d8 

—- = lim 
ds a 



■ exp 



s2 



(56) 



Taking into account the finite resolution of the simulation grid, a finite value of as 
is necessary despite the desired limit ^ 0. In the calculation we use as = Ap 
and choose the resolution of the momentum jump grid as As = 0.1 Ap in order 
to resolve the structure of d6/ds in the s integration. A further reduction of as 
(and an according refinement of the momentum jump grid) does not increase the 
quality of the results. In addition, the maximum allowed time step has to be severely 
reduced to prevent numerical instabilities due to the increased uj{s^ qr) values in the 
vicinity of 5 = 0. For the relevant range of s and qr the function u;{s^qr) is shown 
in Fig. [sjfor Vi{q) to V4{q). Using the same phase space grid as for the first order 
approximation, the time step At = t — to = 0.04ii^ fulfills the stability requirements 
for the chosen value of as = Ap. 

Tomographic approach The crucial point of the tomographic approach is a suit- 



able sampling of the potential landscape entering (29). A straight forward imple- 



mentation of (30) suggests the consideration of the whole coordinate axis for each 



time step and each trajectory (X, /i, v). Then depositing each varied Gaussian onto 
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Fig. 8. Weighting function a;(s,gr) in ( |22| for the potentials Vi{q) — V4{q) (left to right). The 
shown functions u;(s,qr) are normalized to uUi^UmUe = 300 (450) for the potentials V{i 2}(Q') ( 
^{3,4} (9))- The derivative of the ^-distribution is implemented according to ( |56| > with as = Ap = 
0.045np. Note that the oscillation period in s direction is influenced by the range of considered qr 
and q' in the integration ( |22| ). Here we used \q' /u£\, \Qt/u£\ < 30,30,30, 10. 



the coordinate axis to obtain such an implementation would be closely 

related to the concept of Sect. \1.2\ The main difference between both methods is 
the more complicated propagation of the individual trajectories in the tomographic 
approach. Using a local harmonic instead of a linear expansion of the potential and 
including a diffusive term allows this method to use a larger time step. But the 
computational overhead caused by the complexity of the calculation of each time 
step clearly outweighs this profit. 

Exploiting the major advantage of the tomogram - its posit ivity - calls for 
choosing the used coordinates by a MC procedure instead. A direct sampling of the 
coordinates according to the current probability density \ip{q^t)\'^^ however, fails 
completely to reproduce the exact results. Instead of the splitting, only a diffu- 
sive broadening of the initial wave packet is observed for Vi{q). Apparently no 
trajectories overcome the barrier as in (29) only the repulsive force from the po- 
tential but not the actual momentum distribution is taken into account. To include 
the interplay between potential and momentum, for the data presented in Figs. [5] 
and [6] we evaluate the potential at the positions of simultaneously propagated aux- 
iliary trajectories. Starting from a phase space point (q^p), sampled from the initial 
(quantum) state, they are classically propagated in time. Note that the auxiliary 
trajectories determine the potential entering in (30), but the evolved tomogram 
exerts no back action on them. Thus (X, /i, i/), and correspondingly the center of 
the varied Gaussian for the deposition of may significantly deviate from 

the auxiliary trajectory, especially if the potential is evaluated in a region of neg- 
ative curvature. The diverging signatures around t/ut = 5,35 for Vi{q) are due to 
auxiliary trajectories with energies of almost exactly Vq. Those stay in the negative 
curvature region in the vicinity of the potential maximum for a long time. There 
the evolution of (X, /i, u) is governed by hyperbolic functions in ([29|. 

Apart from this, the pronounced smoothing of the results in Fig. [5] is striking. 
The broadening during the time evolution results from the dependency of the width 
of the deposited Gaussians on (X, /i, i/), spoiling a good resolution. On the other 
hand, approximate results are accessible with much less {Np = 12000) trajectories 
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than for the other discussed semiclassical methods. Controhing the steadily increas- 
ing width of the Gaussians requires some kind of restarting procedure in which 
the tomogram is resampled by Gaussians of unit width after a certain time. Per- 
forming such a resamphng at each time grid point, we again approach the concept 
behind Sect. |1.2[ Dechning such a restarting procedure in this work, the given re- 
sults demonstrate the limitation of the tomographic approach when sampling the 
auxiliary trajectories only once. 

3. Conclusion 

In this work, we have compared different semiclassical approaches to quantum me- 
chanics regarding their numerical implementation and efficiency. Focusing on the 
time evolution of a wave packet in one-dimensional quantum structures, we stud- 
ied tunneling, interference and nonlinearity effects. Results were obtained for the 
probability density and various expectation values and contrasted against the exact 
quantum mechanical solution, calculated by means of a Chebyshev expansion tech- 
nique. The Chebyshev method is fast, numerically stable and therefore perfectly 
suited to resolve the full dynamics of a quantum system. Accessible system sizes 
are much larger than the ones that can be reached by other direct solution schemes 
of the time-dependent Schrodinger equation (e.g., using the Crank-Nicholson algo- 
rithm or full diagonalization) . 

A brute force implementation of the Feynman path integral can be performed by 
adapting a linearized semiclassical propagator method, where the inclusion of 'all 
possible paths' is traced back to the set of possible initial conditions on a discrete 
coordinate and momentum grid. Having in mind that within this approach the 
computation time scales as N'^Nf^ where N (Nt) is the number of space (time) grid 
points, the applicability to more complex systems is obviously limited. Instead, the 
implementation should be considered as a 'proof of principle' that quantum effects 
are accounted for correctly if one takes into account the complete superposition 
of the complex weighted trajectories within a local linear approximation of the 
potential. Implementations going beyond this linear potential approximation require 
a full inclusion of the monodromy matrix. If one, along this line, correctly takes 
into account any phase jumps at the focal points, the time step may be increased 
without loss of accuracy. On the other hand, the neglect of some trajectories by 
the MC sampling procedure of the initial conditions leads to a systematic loss of 
accuracy. 

Adopting a probabilistic point of view, the Wigner representation of quantum 
mechanics offers an alternative approach to quantum dynamics. In the Wigner- 
Moyal scheme the Wigner function is propagated in time according to an equation 
of motion, being equivalent to the von Neumann equation. We transform this equa- 
tion of motion into an integral equation which can be solved in terms of an iteration 
series. Then, to leading (first) order, classical trajectories are propagated in time. 
Thereby their initial conditions are sampled from the initial Wigner function. Here, 
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the low computational costs outweigh the loss of some aspects of quantum dy- 
namics. Trying to improve the quality of the approximation by including the next 
(second) order term of the iteration series, we are faced with a tremendous increase 
of computation time. This is due to the necessity of considering a large number of 
trajectories, a high resolution of the phase space grid and a correspondingly small 
time step, in order to avoid the amplification of numerical fluctuations. Thus, in 
order to get the exact quantum mechanical results, we need an even higher numer- 
ical effort than for the linearized semiclassical propagator method. Contrasting the 
Wigner function results of both orders, the gain in accuracy for the second order 
scheme is only moderate such that the additional computational overhead seems 
not justified. Hence, the complexity of the implementation and the ill posed nu- 
merics impede the application of the higher order Wigner- Moyal approach to the 
description of more complex systems. 

Finally, the tomographic representation of quantum mechanics aims at describ- 
ing quantum dynamics in terms of a positive semidefinite function. The quantum 
tomogram can be interpreted as a set of Radon transformations of the Wigner func- 
tion. Relating its time evolution to a diffusive Markov process, the dynamics of 
the system is governed by a set of stochastic integral equations derived from the 
Kolmogorov equations for the tomogram. In the calculation of the drift and dif- 
fusion coefficients for the stochastic differential equations, we used the harmonic 
propagator deduced from a local, second order approximation of the potential. The 
sampling of the potential landscape is the crucial point of this method and strongly 
influences the quality of the data. Evaluating the potential at the coordinates of 
classically evolving trajectories, we reproduce the quantum results qualitatively but 
not quantitatively. As compared to the other considered semiclassical approaches, 
the quality of the data is poor, especially for potentials with distinct negative cur- 
vatures. We expect a noticeable improvement of the results only if a more efficient 
sampling of the coordinates for the potential evaluation can be found. With respect 
to the computational costs, the tomographic approach is slightly more expensive 
than the Wigner- Moyal approach in first order approximation but much faster than 
the linearized semiclassical propagator method or the Wigner-Moyal approach in 
second order. 

To summarize, although the above analyzed semiclassical methods in principle 
capture all quantum effects, they largely differ in quality and required computational 
costs. If one is interested in a method to include minor quantum corrections on 
top of a classical description, the first order Wigner approach is best suited as it 
provides a reasonable compromise between accuracy and computation efficiency. 
Of course, it is possible to reproduce the complete quantum mechanical solution 
at the expense of a dramatic increase of the computing resources. In this respect, 
the linearized semiclassical propagator method is slightly more efficient than the 
second order Wigner approach. While more quantum effects should be included 
in the tomographic than in in the first order Wigner description, the expectation 
values obtained by the former approach are not as good as expected. These aspects 
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have to be kept in mind when applying the above methods to the time evolution of 
more complex systems. 
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